function mean_error = doppler_linearity_error(satpos,satvel,recpos,recvel,d,tu,filename)
%DOPPLER_LINEARITY_ERROR 计算doppler线性化误差,输入都是ecef下坐标
%   此处显示详细说明
velstation = recvel;
posstation = recpos;
x_r = posstation(1);
y_r = posstation(2);
z_r = posstation(3);

chose_sat_pos = satpos;
chose_sat_vel = satvel;
% x_s = chose_sat_pos(1);
% y_s = chose_sat_pos(2);
% z_s = chose_sat_pos(3);

%计算线性化误差
f_vr = (chose_sat_vel-velstation)*(chose_sat_pos-posstation)'/norm(chose_sat_pos-posstation);

d_error = [];
f_error = [];
n = 1;
%% 地心与测站连线上，线性化误差与初始点到测站的距离d的关系
for i = 0:0.001:1
    pos_initial = i*posstation;
    x_i = pos_initial(1);
    y_i = pos_initial(2);
    z_i = pos_initial(3);
    d_error(n) = norm(posstation-pos_initial);%迭代初始点与测站真实位置的距离
    P_ini = norm(chose_sat_pos-pos_initial);
    delta_v = chose_sat_vel-velstation;
    delta_r = chose_sat_pos - pos_initial;
    df_dr = (delta_r.*(delta_r*delta_v'))/P_ini^3-delta_v./P_ini;
    f_vr_ini = (chose_sat_vel-velstation)*(chose_sat_pos-pos_initial)'/norm(chose_sat_pos-pos_initial);
    df_dx = df_dr(1);
    df_dy = df_dr(2);
    df_dz = df_dr(3);

    p_linear = f_vr_ini+df_dx*(x_r-x_i)+df_dy*(y_r-y_i)+df_dz*(z_r-z_i);
    f_error(n) = abs(f_vr-p_linear)/f_vr;
    n = n+1;
end
% figure(1)
% plot(d_error,f_error);
%% 初始点到测站的距离d一定的情况下，初始点在固定球面上时，线性化误差的变化情况

f_error2 = [];
n = 1;
li = 0.01;
bound = 2;
for beita = 0:li:bound*pi
    for alpha = 0:li:bound*pi
        delta_ini = [d*cos(beita)*cos(alpha),d*cos(beita)*sin(alpha),...
            d*sin(beita)];
        pos_initial = posstation+delta_ini;
        x_i = pos_initial(1);
        y_i = pos_initial(2);
        z_i = pos_initial(3);
        P_ini = norm(chose_sat_pos-pos_initial);
        delta_v = chose_sat_vel-velstation;
        delta_r = chose_sat_pos - pos_initial;
        df_dr = (delta_r.*(delta_r*delta_v'))/P_ini^3-delta_v./P_ini;
        f_vr_ini = (chose_sat_vel-velstation)*(chose_sat_pos-pos_initial)'/norm(chose_sat_pos-pos_initial);
        df_dx = df_dr(1);
        df_dy = df_dr(2);
        df_dz = df_dr(3);
    
        p_linear = f_vr_ini+df_dx*(x_r-x_i)+df_dy*(y_r-y_i)+df_dz*(z_r-z_i);
        % err_surf.beita(n) = beita;
        % err_surf.alpha(n) = alpha;
        err_surf.f_error2(n) = abs(f_vr-p_linear)/f_vr;
        % f_error2(n) = abs(f_vr-p_linear)/f_vr;
        n = n+1;
    end
end
beita = 0:li:bound*pi;
alpha = 0:li:bound*pi;

f_error2 = abs(err_surf.f_error2);
f_error2 = reshape(f_error2,[length(beita),length(beita)]);
mean_error = mean(mean(f_error2));
disp(mean_error);
% figure(2)
if tu == 0
    return;
end
if tu == 1
    figure
    mesh(beita,alpha,f_error2);
    xlabel("elevation (rad)", 'rotation', 20);ylabel("azimuth (rad)", 'rotation', -20);
    zlabel("error", 'rotation', 90);
    title('linearization error');
    plot_beautier;
    if filename==0
        return
    end
    print('-dpng','-r200',filename);
end
if tu == 2
    figure
    gca=pcolor(beita,alpha,f_error2);
    colorbar
    set(gca,'linestyle','none');
    xlabel("elevation (rad)");ylabel("azimuth (rad)");
    title('linearization error');
    plot_beautier;
    if filename==0
        return
    end
    print('-dpng','-r200',filename);
end
% f_error2 = reshape(f_error2,[63,63]);
% pcolor(beita,alpha,f_error2);
% colorbar
end

